Self-assembly mechanism in colloids: perspectives from Statistical Physics 
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Motivated by recent experimental findings in chemical synthesis of colloidal particles, we draw an 
analogy between self-assembly processes occurring in biological systems (e.g. protein folding) and a 
new exciting possibility in the field of material science. We consider a self-assembly process whose 
elementary building blocks are decorated patchy colloids of various types, that spontaneously drive 
the system toward a unique and predetermined targeted macroscopic structure. 

To this aim, we discuss a simple theoretical model - the Kern-Frenkel model - describing a 
fluid of colloidal spherical particles with a pre-defined number and distribution of solvophobic and 
solvophilic regions on their surface. The solvophobic and solvophilic regions are described via a 
short-range square-well and a hard-sphere potentials, respectively. 

Integral equation and perturbation theories are presented to discuss structural and thermody- 
namical properties, with particular emphasis on the computation of the fluid-fluid (or gas-liquid) 
transition in the temperature-density plane. 

The model allows the description of both one and two attractive caps, as a function of the fraction 
of covered attractive surface, thus interpolating between a square- well and a hard-sphere fluid, upon 
changing the coverage. 

By comparison with Monte Carlo simulations, we assess the pros and the cons of both integral 
equation and perturbation theories in the present context of patchy colloids, where the computa- 
tional effort for numerical simulations is rather demanding. 
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I. INTRODUCTION 

Self-assembly is a process in which components spontaneously form ordered aggregates^. There exist several exam- 
ples of this mechanism, as it is ubiquitous at vary different length and energy scales. 

Under appropriate conditions, atoms form crystals that are kept together by strong covalent bonds^. This occurs 
because the crystal is the ground state of the system with a large energy gap compared with any other disordered 
aggregate of the same atoms, and hence the crystal structure is very robust with respect to any possible defect 
formation. The specificity of the interactions and of stoichiometry rules, largely constrains the number of possible 
final structures that can be obtained, so that only a relatively small number of crystal structures are found. 

At the opposite side stand soft matter systems where interactions are much weaker, of the order of thermal energy^. 
This includes colloids, polymers, and micelles that also can form ordered or disordered macroscopic aggregates under 
suitable conditions. In this case, however, bonds are easily broken and reformed in a solvent, thus allowing the system 
to rearrange itself in many different ways with configurations separated one another by small energy barriers. The 
combined effect of weak interactions and low specificities, hence results in the possibility of obtaining many final 
aggregates, typically plagued by many defects and low reproducibility. 

Biological systems make an exception to this rule. For instance, a one-dimensional string of amino acids folds 
reliably and reproducibly into a well defined unique native stated, in striking contrast with a similar collapse tran- 
sition characteristic of a homopolymer, where the globular structure has many local minima with very small energy 
differences. The reason for this very different behavior of the polypeptide chain can be traced back to the specificity 
of the amino acid sequence combined with the steric effects provided by the side chains^. 

A similar philosophy has been recently followed in engineering the formation of macroscopic well defined target 
structure starting from decorated (patchy) colloids having multiple functionalities and hence specific interactions^. 
This has become possible in view of the tremendous improvements that have been achieved in the chemical synthesis 
of colloidal particles 7 -^. 

While standard colloids display the typical behavior discussed above^, with many energetically similar aggregated 
structures, patchy colloids are also driven by van der Waals-like forces to a spontaneous aggregation in a huge variety 
of ways, but with only one ground state that is dictated by the specificity of the anisotropic interactions^. 

A system of patchy colloids, therefore, appears to be one of the most promising route to achieve sufficient control 
of the design and the interactions between the building blocks, and hence mimic biological systems such as protein 
folding^, a very appealing perspective from the technological view point to create very complex ordered structures at 
the micro- and nanoscalei. 



2 



The fact that this is indeed a realistic possibility has been recently showed by a set of experiment s] 10 ' 11 , supported 
by numerical simulation s) 12 ' 13 . Here spherical colloidal particles having two hydrophobic poles separated by a charged 
inner band, were observed to self-assembly into a two-dimensional well defined ordered structure - not the most natural 
one associated with the symmetry of that particular colloid, in a specific region of the phase diagram. This will be 
recalled in Section [TT1 

In the remaining of this paper, we will discuss the theoretical approach to this problem in the framework of a 
specific model of patchy colloids, the Kern-Frenkel modelM(see Section HIT)) , that is the same used to numerically 
reproduce the above experimental finding :) 12 ' 13 . In this paper, we will focus on some theoretical approaches that can 
capture the essential features of the system behavior, in particular the part of the phase diagram associated with the 
fluid-fluid (or gas-liquid) transition. 

Two schemes will be discussed, both hinging on classical concepts in statistical physics. 

The first one is an adaptation to patchy colloids of the usual integral equation theory^—, and is patterned 
after a similar technique applied to molecular fluids 1 ^ - — . This will be discussed in Section IIV A[ The second is a 
thermodynamic perturbation theory, again originally devised for isotropic potentials 2 ^ 3 -, and more recently applied to 
the Kern-Frenkel potentia l 24 ' 25 (see Section llV B[) . Finally, Section|V]will draw some conclusions and open perspectives. 

II. SELF-ASSEMBLY MECHANISMS: FROM PROTEIN TO PATCHY COLLOIDS 

A. Self-assembly in proteins and in molecules 

A paradigmatic example of self-assembly process in biology is given by the folding of a proteirA Proteins are 
biopolymers formed amino acids that are linked together into a polypeptide chain by covalent bonds. In order to 
perform the function that are designed for, each protein must be in its native state, that is in a globular form with a 
specific distribution of atoms in space. The protein then spontaneously folds, following the typical pattern sketched 
in FigLll from an unfolded (extended) initial state (U), to a well defined specific native state (N), where the system 
is rather compact in space, usually passing through one or more transition states (TS), and some local minima. 




N 

FIG. 1: A typical folding process. An unfolded polypeptide chain (U) folds into a well defined native state (N), through a 
sequence of transitions between intermediate transitions states (TS) and some local minima. 

This self-organization of the protein where a unique target final structure (the native state) is eventually reached, 
although the system might pass through a series of intermediate metastable higher energy states, largely differs from 
the corresponding transition in standard polymers. In the polymer case 2 ^ a large number of possible final globular 
states, each having comparable, if not identical, energies, are conversely possible. In both cases, the energy barriers 
are of the order of the thermal energy and hence a global rearrangement is possible, and this explains the "glassy 



3 



structure" resulting in the polymer case. In the protein case, however, the system is driven toward a well defined 
target structure by some specificities present in the sequence of the amino acids, and by the steric hindrance given 
by side chains that drastically reduce the number of allowed configurations. This allows a global rearrangement but 
within a much lower number of local minima, so that the system eventually manages to find its way toward the native 
state within biological times (of the order of milliseconds) . 

Consider now the opposite case of atoms or molecules that assemble into a class of well defined crystal structures^. 
In this case, the driving forces are strong covalent bonds and the final number of possible structures is small because 
those are the only ones allowed by the stoichiometry of the interactions. The combination of the strong interactions 
(that do not allow for any rearrangement due to thermal energy) and specificity (given by chemical constraints) , forces 
the system to self-assembly into a unique final structure compatible with the specificities of all given constraints. The 
number of possible crystal structures is however small and enumerated in crystallographic studies. 



B. The colloidal domain 



Colloids are mesoscopic particles with a diameter typically varying between 10~ 8 to 10~ 4 meters, dispersed in a 
microscopic fluid^. The latter is formed by atomistic particles (with dimensions of the order of 10~ 10 meters) and is 
then often treated as a uniform continuum, so that simple model potentials (and the corresponding phase diagrams) 
that are used for atomistic fluids are also valid for colloidal suspensions, but at different length scales, as sketched in 
FigfJ 



Atomic scale 



Colloidal scale 



Human scale 



Astronomic scale 



-15 -10 

10 - 10 



-9 -6 

10 - 10 



io" 3 -io 3 



(meters) 
>10 



FIG. 2: The colloidal scale of the order of 10 9 (nano) and 10 6 (micro) meters. 



One of the main advantage of colloidal suspensions, with respect to their atomistic counterpart, lays in the exper- 
imental possibility of an almost arbitrary tunability in the interactions, in striking contrast with molecular systems 
where possible interactions are severely constrained by chemical requirements. An additional attractive feature of 
colloidal particles, especially from the theoretical point of view, is their almost perfect sphericity, so that it is possible 
to set up a colloidal system behaving as hard-spheres for all practical purposes^. 

Conversely, the main disadvantage of colloids, not present in atomistic systems, is related to the fact that particles 
may have slightly different sizes, charges, or other chemical characteristic, a feature known as polydispersity^. 

While this constitutes a problem in the general interpretation of experimental results, several improvements in con- 
trolling polydispersity have been achieved over the years, so that a number of systems now exists where polydispersity 
can be reduced to values as low as 10% of the average particle sizes, and the system can be reckoned as monodisperse 
for almost all practical purposes. 

As other soft matter systems, such as polymers and micelles, colloids have interaction strengths of the order of 
thermal energy controlling their organization in clusters. When combined with non specificity in the interactions, this 
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allows for a countless number of possible complex structures (i.e. target structures), that differ from one another for 
very tiny free energy differences, that is the opposite of what one would like to achieve in a mechanism mimicking 
protein folding. 

As we shall see in the next section, this last shortcoming can be now overcome with patchy colloids, where the 
number of possible target structures is again huge but, in view of the specificities introduced by the anisotropic 
potential, they have energy absolute minima well separated by metastable states having significantly higher energies. 
Patchy colloids then appear to combine the best of the two worlds, a feature making these systems very attractive for 
technological applications, as remarked. 

C. Experimental framework, new chemical synthesis, technological applications 

Chemical synthesis of colloidal particles is far from being new^, but it has recently experienced a resurge of interests 
in view of the remarkable improvement in controlling the composition of the particle surface^. 

Today, there exists a well defined set of experimental protocols to produce in a controlled way a wide variety of 
colloidal particles with different anisotropics in shape, chemical composition, and functionality. In particular, it is now 
possible to directly modify the surface pattern of colloids to produce multiple surface functionality, an experimental 
technique known as templating. This possibility is extremely appealing for technological applications in such, for 
instance, crystalline arrays with complex unit cells could be used to construct materials having photonic band gap 
in the visible^. More generally, this would provide the possibility of designing a bottom-up experiment to achieve a 
predefined target structure starting from a suitable set of elementary units, as sketched in Fig|3] 

FIG. 3: A cartoon of the self-assembly process from the elementary units to the final target structure. Different colors and 
shapes represent the different building blocks (three in the present case) involved in the process. 

A glowing example of this possibility is provided by recent experiments by Granick and collaborators^, where using 
triblock Janus particles, to be defined below, they succeeded in obtaining a particular planar crystal structure known 
as Kagome lattice. 

The triblock Janus particles (see Fig. |4] right) are colloidal particles having the surface partionated in three parts: 
two solvophobic caps at the opposite poles of the sphere (the light region of Fig[J| and one solvophilic internal band 
(the darker region). The name Janus refers to the fact that the sum of the areas associated with the solvophobic 
and solvophilic regions are almost identical. The most remarkable feature of this experiment lays in the fact that the 
triblock Janus particle has not the most natural symmetry to obtain a two-dimensional Kagome lattice, that requires 
each particle to have four nearest-neighbors. A more natural symmetry would clearly be a sphere with four solvophobic 
patches. However experimentally this pattern is found to be extremely hard to obtain in a reliable and controlled 
way. The authors of ReffTTI. supported by successive numerical simulation o 12 ' 13 , nonetheless found that the same final 
target structure can be achieved by using a triblock Janus fluid within a specific region in the pressure-density plane. 

In the original Janus particles (Figj4] left) the surface of the sphere is partionated in only two rather than three 
parts, and this has very interesting features on its own right, as we will see in the next section. They can be obtained 
by using a technique known as Pickering emulsion, that is a clever way of trapping colloidal particles at the interface 
of a water-wax emulsion so that they can be conveniently functionalized on one side only. The advantage of this 
technique, compared with traditional templating techniques where colloidal particles laying on a plane are templatcd 
and functionalized on one side, is related to the higher yield in terms of particle production, due to the lack of the 
two-dimensional constraint. 

Janus particles in pure water tend to repel each other, in view of the strong solvophilic forces acting on them. 
However, this repulsion can be screened by adding suitable salt, thus leading to a spontaneous formation of clusters 
of various geometries^. 

Below, we will address a theoretical model what is able to mimic this behavior. 



FIG. 4: The Janus particle (left) and the Janus triblock particle (right). The light and dark regions correspond to the 
solvophobic and solvophilic regions, respectively. In both cases, the sum of the solvophobic surface areas equals the solvophilic 
one. 



III. THE KERN-FRENKEL MODEL OF PATCHY COLLOIDS 



The model we are considering is due to Kern-FrenkeUi and the idea is the following. We consider a fluid of 
spheres where the surface of the sphere is divided in two parts having square-well and hard-spheres character, the 
first mimicking the solvophobic region, the second the solvophilic, within an implicit solvent description. 

The position of each particle in space is identified be a set of vectors r,, with i = 1, . . . , N N being the total number 
of particles, whereas the augular orientation each square-well patch on the surface sphere is identified by unit vectors 
rii. Finally, the direction connecting centers of spheres i and j are characterized by unit vector fy = (r, — rj)/|r,- — r^. 
Figure [5] depicts the situation in the case i = 1 and j = 2. 

Therefore, two spheres of diameter a attract each other via a square-well potential of width (A — l)er and depth e if 
the directions of the patches on each sphere are within a solid angle defined by 6q with their relative distance within 
the range of the attractive well, and repel each other as hard spheres otherwise. As the system is still translational 
invariant, the pair potential depends upon the difference = — Yj, rather than and Yj separately, and has the 
form 



$ (ry,ni,rij) = <£ H S (nj) + <few (Uj) * (n», nj, ?y) 
where r,j = |r,j -|. The first term in Eq.|T|) is the hard-sphere contribution 



(1) 




< r < a 



a < r 



(2) 



while the second term can be factorized into a isotropic square-well tail 




a < r < Act 
Act < r 



(3) 
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FIG. 5: The one-patch Kern-Frenkel potential, where f 12 is the direction joining the two centers. Directions of the patches are 
specified by unit vectors ni and n.2 and the surface coverage is x = 0-8. 



multiplied by an angular dependent factor 

f 1, if n 4 • f y > cos 6»o and -n^ • f y > cos6l 
w^n,,^-^ otherwise 1 > 

The unit vectors nj(wi), (i — 1,2) are defined by the spherical angles a;, = (#,-, ipi) in an arbitrarily oriented coordinate 
frame and fij(fi) is identified by the spherical angle f2 in the same frame. Here (3 = {k^T) , fee is Boltzmann's 
constant, and T the absolute temperature. Reduced units, for temperature T* = fceT/e, pressure P* = /3P/ p, density 
p* = pa 3 , and chemical potential p* = /3/x are routinely used in the description of the thermodynamic. The above 
potential then ensures that, in order to be bonded, two particles have to be within the range of the square-well 
potential with their attractive surfaces properly aligned. If not, they behave as hard spheres. 

The relative ratio between attractive and total surfaces is the coverage x that is related to the semi-angular width 
80 of the patch. This can be obtained a o 14 ' 16 

X = (^(^n^Z^sm 2 ^. (5) 

where we have introduced the angular average 

(...)„ = -L/du,... (6) 

IV. THE TOOLS OF STATISTICAL PHYSICS 

As discussed above, it is crucial for self-assembly process to have a full control of the phase behavior of the system 
to make sure that a particular condensed phase is obtained for a given set of parameters. 
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/ \ Fourier transform A / 7 \ 

c(r) > C{K) 



Closure 



/ N Inverse Fourier transform ~ / 7 \ 

7(r) < j[k) 



OZ equation 



TABLE I: Schematic flow-chart for the solution of the OZ equation for isotropic potentials. 



Statistical physics has developed a number of different techniques able to assess thcrmophysical properties of simple 
and complex fluids. Here we will mention two of them that have proven to be useful in the context of self-assembly 
patchy colloids. 



A. Integral equation theories 



In the case of spherically symmetric potentials, the strategy to infer thcrmophysical properties of a fluid of density 
p has a long and venerable tradition in integral equation theories^. The basic idea is to be able to compute the total 
correlation function h{r) that is directly linked to all structural and thermodynamical properties of the fluid. In order 
to do that, however, the knowledge of the direct correlation function c(r) is also required. The two are related by the 
exact Ornstein-Zernike (OZ) equation 

h(r) = c(r)+p J dr' c (/) h (|r - r'|) (7) 

but an additional equation is required in order to close the system and compute the two functions. The second relation 
also involves the potential 4>{r), and has the general form 

c(r) = exp [-/?<£ (r) + 7 (r) + B (r)] - 1 - 7 (r) (8) 

where 7 (r) = h (r) — c (r). 

Note that h(r) = g(r) — 1, g(r) being the radial distribution function that is usually computed in numerical 
simulations. The crucial point in this matter, is that this second equation ([8]) cannot be computed exactly and an 
approximation (closure) is required to compute the bridge function B(r), that is known only as an infinite power series 
in density whose coefficients cannot be readily calculated. All practical closures then approximate B (r) in some way, 
and this gives rise to thermodynamic inconsistencies that are well-known in this formalism^. The numerical solution 
of the system of non-linear equations ([7]) and ([H]) then proceeds through a self-consistent procedure going back and 
forth from real to momentum space until self-consistency is achieved, as sketched in Table I. Clearly the reason why 
it proves convenient to go into momentum space, is related to the fact that the convolution appearing in OZ equation 
([7]) greatly simplifies in Fourier form 

1 — pc (k) 

h(k) and c(k) being the Fourier transform of h(r) and c(r) respectively. Use of the auxiliary function 7(7*) in place of 
h(r) is also more convenient in practical calculations, as illustrated in Table tf^. 

The case of angular dependent anisotropic potentials is far more complex from the algorithmic point of view, but 
the philosophy behind the methodology is identical. It was devised in the frame of molecular fluids^, and more 
recently adapted to the case of patchy colloid o 16 i 17 . Here we just sketch the idea, referring to Rcfs. [l5l - [T7l for details. 
The whole procedure hinges on a remarkable piece of work carried out by Fred Lado in a series of papers in the 
framework of molecular fluids^—. 

The iterative solution of the angular dependent Ornstein-Zernike (OZ) equation plus an approximate closure equa- 
tion again requires a series of direct and inverse Fourier transforms between real and momentum space involving the 
bridge function B(rij, n^, rij), the direct correlation function c(jij, hi, n^), the pair distribution g(rij, hi, hj), the corre- 
lation functions h{vij, hi, hj) = g{vij,hi,hj) — l, and the auxiliary function 7 (r^ , n^, hj) = h(jij ,hi,hj)—c(rij,hi,hj). 

Angular dependence introduces additional direct and inverse Clcbsch-Gordan transformations^ between the coef- 
ficients of angular expansions in "axial" frames, with z = in direct space or z = k in momentum space, and more 
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TABLE II: Schematic flow-chart for the solution of the OZ equation for the Kern-Frenkel angle-dependent potential. See Section 
IIV Al for a description of the scheme 



general "space" frames with arbitrarily-oriented axes. In the case of the Kern-Frenkel potential of Section Mil with 
only one or two patches, the cylindrical symmetry of the angular dependence entitles the use of the simpler version 
of the procedure for linear molecules^ rather than the full anisotropic version^!. The results is illustrated in Table 
II and is the extension of the isotropic case given in Table I. For i = 1 and j = 2, the procedure in real space starts 
an iteration in the axial frame with the current values of the axial coefficients 7z 1 ; 2 m('''i2) and uses an expansion in 
spherical harmonics to construct 7(^12, wi, u)2)- The bridge function _B(Yi2, wi, L02) is similarly constructed in a given 
closure approximation, and the latter is used to obtain c(ri2, uji, 0J2), its axial coefficients C; 1 ; 2m (ri2) and, through a 
Clebsch-Gordan (CG) transformation, the corresponding space coefficients c(ri2'>hhl)' 

The corresponding part in Fourier space is achieved by Hankel transform, and a backward Clebsch-Gordan trans- 
formation to return to an axial frame in k space, yielding Q 1 i 2m (/c). A parallel sequence of operations in Fourier space 
starts with the OZ equation to get 7/ 1 / 2?n (fc), followed by a forward Clebsch-Gordan transformation, and an inverse 
Hankcl transform, to provide 7(7*12; hhl)- A final backward Clebsch-Gordan transformation, brings a new estimate of 
the original coefficients 7 , 1 ' 2m ( , ri2)- This cycle is repeated until self-consistency between input and output 7; 1 / 2m (ri2) 
is achieved as before. 

In the present framework a particular closure denoted as reference hypernetted-chain (RHNC) has been used for both 
the square-well and the Kern-Frenkel potential, with the bridge function approximated by its hard-sphere functional 
form. It has few shortcomings and several advantages, compared to other possible closures, that have been analyzed 
in detail in past worker— . Its main advantage here is related to its ability of computing the chemical potential 
and pressure without invoking additional approximations besides those included in the closure^—. For a fixed 
temperature T*, one can then compute the pressure of the gas (colloidal poor) phase P* and of the liquid (colloidal 
rich) Pj* phases, and the corresponding chemical potentials //* and //*. The fluid-fluid (gas-liquid) coexistence line 
then follows from a numerical solution of a system of non-linear equations 



p; [T* 



Pa) 
P*g) 



P-i {T ,Pi) 



(10) 
(11) 



whose solutions are the gas p* and liquid p* reduced densities associated with the coexistence lines. 

In case of square- well with a width amplitude of A = 1.5, the resulting phase diagram in the temperature-density 
plane is reported in Fig. [S] . 

The results from reference hypernetted-chain integral equation theory (IE) are contrasted with those stemming 
from Gibbs ensemble Monte Carlo simulations by Vega et a£22 showing good agreement in both branches, far from 
the critical region. The difficulties in assessing the critical region, is a well-known feature of all integral equation 
theories, and it is related to the aforementioned thcrmodynamical inconsistencies related to the closure. Apart from 
this, the performance of integral equation theory is rather satisfactory, with the agreement with numerical simulations 
not only qualitative but also quantitative. This suggests that the the theory is able to capture the essential features 
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FIG. 6: The computed phase diagram of a square-well fluid with A = 1.5. Results from integral equation theory (IE) are 
contrasted with Monte Carlo simulations (MC) of Ref. [3^ . 



of the phase behavior. It is also worth recalling that the fluid-fluid transition is known to become metastablc against 
fluid-solid transition below A = 1.25, that is for sufficiently short-range attractions^ 3 .. 

The case of patchy colloids is tackled using the same methodolog y 16 > 17 . Again the reference hypernetted-chain 
integral equation theory (IE) is contrasted with MC numerical simulations by Sciortino et a l 34 ^ 35 in Fig. [7J As 
representative value, a coverage of x = 0.8 with A = 1.5 was selected, so this value is intermediate between the 
isotropic square-well discussed above, and the Janus limit corresponding to half coverage (% = 0.5) that is known to 
display an anomalous phase diagra m 34 i 35 . Unlike the isotropic potential, however, patchy systems allow to have the 
same attractive part distributed in different ways on the surface, and this is included in the two cases also reported 
in FigH 

The first case is when all attractive parts are concentrated into a single patch (IP), while the second case corresponds 
to an attraction distributed into two circular patches at the opposite poles of the sphere (2P). Note that, at half 
coverage \ = 0.5, the two cases correspond to the Janus and the triblock Janus particle, respectively. The phase 
diagram of the Janus case turns out to be particularly interesting, as it is ruled by a competing mechanism with the 
tendency to form a condensed phase on the one hand, and a clusterization process tending to form isolated and non- 
interacting micelles (and vesicles) inhibiting the condensation process, on the other hand. The result is an anomalous 
phase diagram with reentrant vaporization 3 ^. 

Two features of Fig. [7] are noteworthy. First, the numerical simulations clearly indicate that the coexistence 
line of the two patches case always lay above the single patch counterpart at equal coverage. This means that, on 
cooling from the high temperature single phase region, the two patches system tends to phase separate before the 
corresponding one patch case, as it was to be expected in view of the higher valence in the 2P case. Second, the 
consistency between integral equation theory and numerical simulations appears to be significantly better in the 2P 
as compared with the IP counterpart. Again, this is not surprising in view of the fact that the one-patch potential is 
much more asymmetric and hence harder to be described by an isotropic approximation of the bridge function, such 
as that used in the present approach. 

For sufficiently low coverages, one might expect the fluid-fluid transition to become metastable against other tran- 
sitions, thus providing the angular analog of the metastability occurring for short-range attractive square-well inter- 
actions discussed previously. This is indeed the case. For the one-patch potential, extensive simulations 3 ^ 3 ^ indicate 
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FIG. 7: The computed phase diagram of a Kern-Frenkel fluid with A = 1.5 and surface coverage y = 0.8. Results from integral 
equation theory (IE) are contrasted with Monte Carlo simulations (MC) obtained in Ref. [l6| and[l7] both in the case of a single 
(IP) and double (2P) patches. The solid line depicts the results from Barker-Henderson perturbation theory (PT) as given in 
Ref 25. 



that a micellization process prevents the condensation process altogether. In the two-patches case, on the other hand, 
even at low coverages the double valence on each particle allows the formation of a network, so that the condensation 
process can take place. However, a crystalization was shown to occur for coverages below \ ~ 0-3, so that the system 
forms a crystal before being able to condensate^ 7 -. 



B. Barker-Henderson perturbation theory 

The idea behind perturbation theory is that this second term in Eq.([T]) can be treated as a perturbation of the 
first, that is of the hard-sphere term. In the present context, this means that we are starting from a fluid of hard- 
spheres and perturbatively decorating the hard-sphere surfaces with a attractive patches. We now discuss how a 
classic perturbation theory due to Barker and Henderso n 23 ' 36 ! 37 can be adapted to the Kern-Frenkcl model to infer 
thermodynamical propertie d 24 ! 25 . 

In the case of isotropic potential, the original method is due to Zwanzig 2 ^, but the most useful version was given 
by Barker and Henderson few years later—. The method was recently extended to the Kern-Frenkel potentia l 24 ' 25 , 
and it is based on a second-order cumulant expansion of the excess free energy 

In order to illustrate the method, let us derive the first order term. Let us work in the canonical ensemble, although 
it proves much more convenient to work in the the grand-canonical ensemble, as it correctly takes care of the finite 
size contributions^ 6 -. 

Assume the total potential U to have the general following form 

U 7 (1,...,N) = Uo(l,...,N) + -,U I (l,...,N) (12) 

= J2®-y = $ ° (y) + $i (v) ' 

i<3 i<j i<j 
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where Uo(l, . . . , N) = J^i j ®o(ij) is the unperturbed part and Ui(l, . . . , N) — J^i j is the perturbation part. 

Here < 7 < 1 is used as perturbative parameter. Note that when each coordinate i includes both the coordinate 
and patch orientation n^, so that i = (r^ , n^) , the expression is also valid for the Kern-Frenkel mode l 24 ' 25 . For simple 
fluids, instead i = rj only. 

The potential represented in Eq.([T]) can be clearly cast in this form upon combining the two parts of the second 
term appearing in Eq.([T]) 



(Tij,B4,hj) 

Introducing the following short-hand notation 



5S W [Tij 



N 



(r 4J )*(n 4 ,n J; fy) 



for the integration over all particle coordinates, the partition function 



1 



MA T N 



1 



1,...,N 



N\A 3 T N 7 



(13) 



(14) 



(15) 



(here At is the de Broglie thermal wavelength, and Z 1 is the configurational partition function), can then be used to 
obtain an expansion of the Hclmholtz free energy 3 ^. 



F-. 



7=0 



1 2 

2! 7 



d 2 F 1 
dj 2 



(16) 



7=0 



that is valid for arbitrary 7. 

Taking the derivative of F 1 one has, using Eq. ([T3"|) 



d_ 

d-. 



\P*"y] = \J |^[-^7(12)] Py (12) 



where 



Pi (12) 



N\ 



(N-2)\Z J J 3 _ 



(17) 



(18) 



When 7 = 1, this yields the free energy correct to first order. 

An explicit expression of the second order term is far more involved and can be found in Refl25l. In the particular 
case of the Kern-Frenkel potential, one then obtain o 24 ' 25 the excess free energy per particle with respect to the 
corresponding hard-sphere free energy Frs/N to have the reduced form 



/3(F-F HS ) 
N 



N N 



with 



and 



N 



^ = drr2 ^w (r) go (r) {p*(12)) uii 

677 / dr] 



o 3 \dP* 



drr 2 g (r) 2 w (r) [/?* (12)] 



T Ja 



(19) 



(20) 



(21) 



Here Pq = f3Po/p is the reduced pressure of the HS reference system, and go(r) the corresponding radial distribution 
function. We have also introduced the packing fraction rj = npa 3 /6. 

Once that the free energy is known, pressure and chemical potential can be obtained from the following standard 
thermodynamical relations 2 ^ 



(3P 
P 



d f(3F 



V dr] V iV 
r d ( pF 



(22) 
(23) 
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so that the same procedure given in Eqs. (fT0|) and (JTTJ) can then be applied to infer the phase diagram in the 
temperature-density plane. 

In spite of the fact that this is a high temperature expansion - terms flFi/N and PF2/N are proportional to 1/T* and 
(1/T*) 2 respectively, this scheme is found to provide a remarkably good description of the fluid-fluid transition. This 
is displayed again in FigJTl where the fluid-fluid coexistence curves computed using Barker-Henderson perturbation 
theory (PT) again in the case A = 1.5 and x — 0-8 (that is the same case as before) are reported and compared with 
previous results from integral equation theory (IE) and Montecarlo simulations (MC). 

Given the second-order truncation in the high temperature expansion, the performance of the perturbation theory 
is rather remarkable. It is worth stressing, however, that this was also shown to be the case even in the fully isotropic 
square-well potential^. An additional strong point of this approach, presented in Refl25l. is that even the fluid- 
solid coexistence phase is found to be well reproduced by the Barker-Henderson method. In this latter respect, an 
interesting additional improvement could stem from a methodology akin to that presented in Refl39l 

V. CONCLUSIONS AND OPEN PERSPECTIVES 

In the present paper, we have discussed the phase behavior of a particular fluid of patchy colloids and its potential 
applications to self-assembly processes at micro- and nanoscale, of biological inspiration. 

Interactions in patchy colloids were described by the Kern-Frenkel model where circular attractive square-well 
patches decorate the surfaces of hard-spheres. This model then smoothly interpolates between two well-known and 
studied isotropic fluids - a hard-sphere fluid (coverage \ = 0) and a square- well fluid (coverage x = 1.0). It is 
sufficiently rich to display an interesting behavior, and sufficiently simple to be amenable to a detailed theoretical 
treatment. 

Wc have discussed the crucial importance of having a complete control of the phase diagram of the system, including 
the exact location of the fluid-fluid (gas-liquid) and fluid-solid coexistence lines. Compare to standard isotropic 
colloids, the phase diagram of patchy colloids may be far more complex, and this was hinted in the simple case of 
the Kern-Frenkel potential, where the attractive part covers only a fraction of the total surface, that can also be 
distributed among different parts of the surface. This was specifically addressed when the entire square-well part was 
condensed into a single region (one- patch) or partioned in two parts at the opposite poles of a sphere (two-patches). 

In order to tackle the calculation of the fluid-fluid transition, we have briefly outlined two strategies based on two 
well-known statistical physics tools, namely integral equation theories and perturbation theory. 

In the case of integral equations, the reference hypernettcd chain (RHNC) has been used to compute the phase 
diagram in the temperature-density plane, with the fluid-fluid coexistence lines following very closely those from 
numerical Monte Carlo simulations both in the single and double patches case for a given width of the square-well 
and a given surface coverage. Integral equation theory has two main strengths. First, they can be formulated in 
both the single and double patches case with only minor modifications; second it is quantitatively predictive at a very 
limited computational effort as compared to the corresponding Monte Carlo simulations that instead require a major 
computational effort. The main drawback of the method is related to its internal thcrmodynamical inconsistencies 
preventing the access to the critical region, that has then to be extrapolated in some way. 

Thermodynamic perturbation theory, on the other hand, does not suffer from this drawback, as the full critical 
region, including the critical point, can be evaluated. It has however a different disadvantage that we have tried 
to underline. Our formulation builds upon a classic second-order theory due to Barker and Henderson for simple 
fluids, and we have discussed how this can be adapted to the Kern-Frenkel potential in the case of a single patch. 
Again by means of a direct comparison with Monte Carlo simulations on the same model, we have assessed the 
surprising potentialities of this technique that go beyond expectation, and its good performance in its predictions 
makes it a method competitive with integral equation theory. In this case, the main drawback hinges on the fact 
that, in its present formulation, the exact location of the patches position is irrelevant, and therefore cannot be used 
to discriminate between the single or double patches scenario. 

While in principle other similar theoretical methods such as the Weeks-Chandler- Anderson (WCA) perturbation 
theory^Sr— or the reference interaction site model (RISM) 29 ' 43 could be possible, in practise those described in the 
present work and their possible improvements^, are expected to be the most efficient ones, in view of the particular 
all-or-none angular dependence of the attractive part. 

The topic of patchy colloids treated in the present paper is a very active research topic, and many different 
approaches have been envisaged to tackle the analysis of the phase behavior in details. Although an exhaustive list 
of them is clearly beyond the aim of the present work - it can be found in Ref. l44l . we would still like to mention 
some very recent results as they are relevant for the point of the present discussion. Fantoni et al^ used a cluster 
theory to rationalize micelles formation in the Janus fluid. Reinhardt et al^- mapped the effect of the patchy colloid 
anisotropy into a binary mixture of isotropic colloids to explain the Janus anomalous phase behavior observed in 
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Refill A different model with similar aims has also been studied by Tavares et al^. In this model, earlier proposed 
by Bianchi et al^, attractive patches are replaced by sticky spots distributed over hard-spheres surfaces. Here, sticky 
spots act as a specific binding interactions so that the ground state has all possible favorable contacts saturated, with 
hence an energy significantly lower with respect to the other metastable states. However, the point-like nature of the 
attractive spots makes a local rearrangement rather difficult, and the system has then the tendency to be trapped 
into a metastable state, in spite of the large energy gap, a drawback not present in the Kern-Frenkel potential, where 
the number, the extension, and the location of the patches can be freely and independently tuned to create the most 
favorable conditions. 
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